Ground-state energy of the electron liquid in ultrathin wires 
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The ground-state energy and the density correlation function of the electron liquid in a thin one- 
dimensional wire are computed. The calculation is based on an approximate mapping of the problem 
with a realistic Coulomb interaction law onto exactly solvable models of mathematical physics. This 
approach becomes asymptotically exact in the limit of small wire radius but remains numerically 
accurate even for modestly thin wires. 



Recently much attention has been devoted a class of 
one-dimensional (ID) conductors that can be termed 
ultrathin wires. Examples of such systems include 
single-wall carbon nanotubes (CN) semiconductor 
nanowires 0, and conducting molecules Q. Semicon- 
ducting ultrathin wires are especially interesting because 
their electron concentration n can be varied by the field- 
effect, which can be used for creating miniature electronic 
devices 0, 0| . In such applications a fundamental role is 
played by the concentration-dependence of the ground- 
state energy density s(n). This function determines the 
electrostatic screening and affects the capacitive coupling 
of the electron liquid to external voltage sources. It is 
also a core input of the density- functional theory (DFT), 
which is the basis of today's electronic structure calcu- 
lations. Since the 3D Fermi-liquid theory does not hold 
in ID, it is unclear whether the usual DFT optimized 
for three-dimensional (3D) systems is adequate for ultra- 
thin wires. The Luttinger liquid (LL) theory which is 
called upon to replace the Fermi-liquid theory, makes no 
predictions for the short-range physics that determines 
e(n). Therefore, the calculation of the ground-state en- 
ergy of ID wires with realistic Coulomb interactions has 
remained an open problem. The primary difficulty is 
the computation of the correlation energy e cor , which 
is determined by the shape and size of the exchange- 
correlation hole (XCH), i.e., the reduction in probability 
of any two electrons closely approaching each other. Be- 
low we propose a theory that calculates these quantities. 

Model. — Our calculation is done for an A-component 
electron gas, iV being the combined spin- valley degener- 
acy of the electron spectrum. For example, A — 4 in 
CN . The aforementioned XCH is the term that refers 
to the negative dip of the two-body cluster function h(x) 
around x = 0. Here h(x) — (Mn)^ 1 J2i^j(^( x i ~ x j ~ 
x)) — 1, M is the number of electrons, and their 
coordinates. Larger \h(0)\ imply stronger correlations. 
Since /_ h(x)dx = —1/n, the XCH has a characteris- 
tic width Z» ~ l/n|/i(Q)|. For example, in the free Fermi 
gas h(x) = - A sin 2 (nx/N)/n 2 x 2 , so that |/i(0)| = 1/N 
and Z» is equal to N/n, the average distance among elec- 
trons of same species or, as we call it, same isospin. Our 
goal is to compute h(x) for an interacting system. Once 
h(x) is known, e(n) can be obtained straightforwardly, 



see below. We model the interactions by the potential 
U(x) = e 2 1 k\/ x 2 + R 2 , which accounts for smoothing of 
Coulomb repulsion at distances of the order of the wire 
radius R. The wire is considered ultrathin if the parame- 
ter C = ln(as/i?) is large, where as = h 2 K,/me 2 , to, and 
k arc the effective Bohr radius, electron mass, and dielec- 
tric constant, respectively. On general grounds, we may 
expect that at low densities, n -C 1/ag, electrons should 
form a ID Wigner "crystal" Q with h{x) sharply peaked 
at integer multiples of a — 1/n. At n 3> I/as where 
electrons have a large kinetic energy, h{x) should remain 
appreciable down to i < a. Below we refine and flesh 
out this qualitative picture by quantitative calculations. 

Crucial for our approach is the fact that to the leading 
order in 1 / C the problem in hand and the problem with 
the contact interaction, U(x) = (h 2 c/m)S(x) give the 
same short-range behavior of the correlation functions, 
including the XCH. Here c is given by 



c= (2/a B )Hh/R). 



(1) 



This remarkable mapping between the two interaction 
laws holds only in the liquid state, n> 1/as- The reason 
for it becomes clear if one carefully separates the effects 
of the sharp maximum ( "core" ) of the Coulomb potential 
U(x) at x = from those of its 1/x tails. As was shown in 
our earlier work 9], the condition n» 1/ob guarantees 
that the Coulomb tails have negligible effects on h(x) up 
to exponentially large distances, \n{x/a) ~ l/r s , where 
r s = a/2as <C 1. Since r s plays the role of the dimen- 
sionless coupling constant, this agrees with the conven- 
tional wisdom. On the other hand, the electron scatter- 
ing caused by the short-range core of U(x) is enhanced 
by the large logarithm C. Therefore, the Coulomb poten- 
tial acts as a sum of a strong short-range core and weak 
tails, and so can be mapped onto a suitable (S-function. 

The rest of the paper is organized as follows. We be- 
gin by studying certain limiting cases, which verify the 
correctness of our choice Q of the coefficient c. We then 
explain how our theory can be used to calculate e(n) 
at all r s <C 1. We proceed to the study of the large- 
r s Wigner crystal where the mapping onto the contact 
interaction model is no longer valid. We show that the 
exact asymptotics of e(n) in the r s 3> 1 limit can nev- 
ertheless be derived while at r s ~ 1 a simple variational 
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approximation can be used. We also present a numeri- 
cal scheme that unifies all the asympotical formulas we 
obtain. It yields a seemless interpolation over the entire 
range of n even for C ~ 2 — 3. We take it as evidence that 
our theory remains numerical accurate even for modestly 
thin wires, which may stimulate its use in practical DFT 
calculations. (Achieving large £ is feasible [tj but tech- 
nically difficult). 

Definitions. — We do the usual subtration of the Hartree 
term in the definition of the energy density, e(n) = 
L~ 1 [(H) -U{{))n 2 /2] , where L is the length of the wire, H 
is the Hamiltonian, and tilde denotes the Fourier trans- 
form. We further define the correlation energy density 
e cor as the difference between e and the sum of the ki- 
netic £q and the exchange e x energies of the Fermi gas, 
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where At — §— 7— hi7r w —0.222, 7 is Euler constant [Io|. 

The relations among £ cor (n), h(x), and the dielectric 
function e(q,u) are (see, e.g., Ref. 0, Sees. 5.4 and 5.6): 



n3 f dr^e int (n,r s , N) ^ 



£int = n 2 / dxU(x)h(x), 
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RPA regime. — The validity of our mapping between the 
Coulomb and the contact interactions can be verified by 
an independent method if the limit of large N (actually, 
large iV 2 ) is taken. We discuss it because it is not only 
an instructive example but also the case relevant for CN, 
where N 2 = 16. For large N, e(q,ui) is dominated by 
the random-phase approximation (RPA) |7j, which sums 
order-by-order the diagrams with the largest number of 
fermion loops. For q ^> fop = im/N the result is 



e(q,u) = 1 + 



2nE{q)U{q) 



E(q) = 



h 2 q 2 



(0) 



E 2 (q) - (hu -MO) 2 ' 2m 

Combined with Eq. JSJ, it entails that at C <C nan <C 
N 2 C (the RPA regime), the XCH has the depth \h(0)\ ~ 
(nnh) -1 and a characteristic width 

I* = y/a B /2nln{l*/R). (7) 

The XCH is much deeper than in the Fermi gas, \h(0)\ 3> 
1/N, and so the correlations are strong; yet \h(0)\ <C 1, 
so that the RPA is still reliable. From Eq. JJJ we find, 
to the leading order in 1/N, 
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Repeating the same calculation for the contact interac- 
tion with c given by Eq. JIJ, we obtain exactly the same 
result. To track down how this comes about it is conve- 
nient to do the integration in Eq. J3J in the g-space. The 
interaction potential enters through its Fourier transform 
U(q) ~ [2e 2 / k) \nx{l/qR), which is a slow function of q. 
The integral is dominated by q ~ l/l*, and so to the lead- 
ing order in U{q) can be replaced by the U{l/U), 
i.e., U{x) — * (ti 2 c/m)5(x), as we claimed above. 

The RPA eventually breaks down at small n, where 
it predicts h(Q) to drop below the strict lower bound of 
— 1 required by the nonnegativity of the electron density. 
This places the lower boundary of the RPA regime at 
n ~ C/iraB- What happens at lower n is discussed next. 
CTG regime. — The case of n <C C/as has in fact al- 
ready been studied in Ref . |IJ . We showed that at such n 
electrons should form a correlated state of the Coulomb 
Tonks Gas (CTG). The CTG can be defined as the state 
where on short lengthscales electrons behave as impene- 
trable but otherwise free. It owes its name to a certain 
similararity it enjoys with the Tonks-Girardeau gas of ID 
cold atoms ^lj- It i s worth mentioning that the long- 
distance behavior in the RPA, CTG, and Wigner crystal 
regimes is universally the same and is described by the 
LL theory. In the limit R — > +0, i.e., c, C — > 00, the 
ground-state wavefunction ^ factorizes into the isospin 
part $ and the orbital part (the remainder) 0: 



(9) 



where Ql through QM are the indices in the spatially 
ordered list of the electron coordinates < xq\ < 
. . . < xqm < L (periodic boundary conditions are as- 
sumed), (— 1)^ is the parity of the corresponding per- 
mutation, and A = 1 for now. For N — 2, $ coincides 
with the ground-state of a spin-1/2 Heisenberg chain; 
for N > 2, see Ref. We will not discuss the func- 
tion W here because it has negligible effect on h(x) for 
x <C aexp(l/r s ) 9]. Once W is set to zero, $ becomes 
the ground state of the contact-interaction problem at 
c = 00 (the gas of impenetrable fermions) \l4 . Il5j . This 
is another explicit demonstration of our mapping, this 
time in the C — > 00 limit. Note that the XCH has the 
largest possible depth of unity and the width h = a. 

For a finite R, ^ remains the correct approximation 
to the ground state to the leading order in £ _1 ,r s <C 
1. We use ^ (with W — 0) as a trial state to evaluate 
e(n). Independent of the form of <I>, the result is given 
by Eq. with N = 1 (see also Ref. 0), 
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which agrees with e(n) for the contact-interaction prob- 
lem to the order 1/C [Eq. I|15fl ]. validating our mapping 
once again. 



Bethe ansatz. — The most remarkable consequence of 
the mapping between the Coulomb and the contact- 
interaction models is that a unified treatment of all 
r s <C 1 regimes is possible. This is due to the fact fact 
that the latter model is solvable by the Bethe ansatz Jl2| . 
The exact energy density at any given n is given by [15j 
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where p{k) is the solution of the integral equation 
Q 

i- = p(k)- f dk'G(k-k')p(k'), 
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Here ^(z) is the digamma function [lOj and Q — Q(n) is 
fixed by the constraint n = J_Qdkp(k). Two analytical 
asymptotics of the solution can be obtained [lil flif 
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c<n<ciV 2 ,(14) 
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in agreement with Eqs. JSJ and (|10fl . From the theory 
point of view, Eqs. @J, Q, |HIJ), and (CH) — (H3I> solve the 
problem of computing e(n) at all r s <C 1. A practical 
algorithm for finding the solution is given shortly below. 
Wigncr crystal. — At very low densities, r 8 3> 1, the 
mapping onto the contact interaction problem is how- 
ever invalid. The tails of the Coulomb barriers that sep- 
arate nearby electrons are strong enough to keep them 
at almost equidistant positions (although the long-range 
order is eventually destroyed by fluctuations). Accord- 
ing to the standard strong-coupling perturbation theory, 
the ground-state energy in this regime is equal to the 
Madelung sum plus the zero-point phonon energy, 



£ = — n 2 [\n(Rn) ~ A w ] + C ph — n 5 / 2 ^ 2 



(16) 



where A w = In 2-7 ps 0.116 and C ph ~ 1.018. As for the 
cluster function h(x), it can be obtained by interpolating 
between the collective phonon-like correlations at x > a 
and two-body correlations at x <C a (see, e.g., Refill 7T. 
Variational and numerical interpolation — Until now we 
expanded on a formalism that gives results both for e(n) 
and for h(x) that are rigorously correct to the leading 
order in a suitable small parameter, either 1/C or l/r s . 
These results, e.g., the functional form of e(n) in various 
regimes [Fig. ^a)] have an academic or methodological 
interest. In the remainder of the paper we shift the focus 
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FIG. 1: (a) Qualitative behavior of e(n). (b) e in units 
of e 2 / K,a 2 B for as/R = 10,15,20 (top to bottom), evaluated 
numerically, (c) The low-density part of the same plot; solid 
lines are from the variational method, the dashed lines - 
from the Bethe ansatz. 



to a more pragmatic goal. We wish to find a computa- 
tional scheme that gives an accurate numerical approxi- 
mation to the same quantities when neigher C nor r s are 
truly large. We achieve this by combining a variational 
method with a numerical interpolation. Some results are 
shown in Fig. ^b) and (c). The concrete interpolation 
scheme used in generating these plots is as follows. For 
r s < 2, e(n) is calculated by numerically solving Eqs. JQ), 
(|TTll-l|lSI) with U = nexp(A T )/Q(ri, c). Note that func- 
tion Q(n) has the following limiting forms: Q ~ im for 
n < c, Q ~ l^fnc for c < n < cN 2 , and Q = nn/N 
for n ^> cN 2 . This entails that our choice of /* is ex- 
act at small and large n, and is adequate everywhere in 
between, see Eqs. ©,0, and (fTU|) . 

All what remains is to handle r s > 2 regime where 
the crossover between the CTG and the Wigner crys- 
tal occurs. Our solution is to treat A in Eq. @ as a 
variational parameter. This ensures a smooth transition 
from the CTG (A = 1) to the Wigner crystal (A > 1), 
provides a strict upper bound on e(n), and can be done 
semi-analytically. Indeed, the energy density of the state 
"J (A) is the sum of the kinetic and the potential e v °^ 



terms. By virtue of a formula similar to Eq. J3Jl, 
can be computed differentiating the known energy den- 
sity £cs(A) in the Calogero-Sutherland model [18| with 
respect to its coupling constant c = A(A — 1), 
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To get e"° t r , we calculate it at A 



i 1,2, and 00 using the 
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exact cluster functions h(x) [19j and interpolate between 
the obtained four values by a cubic polynomial in A -1 , 
e™r = ( e 2 / K )n 2 [ln(i?n)-A w -aiA- 1 -a 2 A- 2 -a 3 A- 3 ]. 
For example, in the C — > oo limit we find a\ ~ —0.3173, 
a 2 ps -0.02363, and a 3 pa 0.003048. Smallness of a 2 and 
a3 implies a high numerical accuracy of this polynomial 
fit. Minimizing e™^ + e™ t r with respect to A (numeri- 
cally), we get e(n). The quality of our variational method 
can be judged by how well it compares with Eq. (|16(l in 
the r s ^> 1 limit. It is easy to see that the functional 
form of e(n) is reproduced correctly, but the coefficient 
in front of the phonon term is approximately 1.022, i.e. 
higher than C p h by mere 0.4%. The results of this proce- 
dure, implemented for several values of cib/R, are plotted 
in Fig.[I](b) and (c). The curves produced by the Bethe 
ansatz and the variational method match virtually seam- 
lessly. Thus, the proposed scheme gives a theoretically 
well-founded and numerically accurate DFT needed in 
applications, some of which are discussed next. 
Implications. — The main physical omission of our theo- 
retical model is the screening of Coulomb interactions by 
other ID subbands that may be present in a wire. Such a 
screening is averted if k exceeds a certain threshold Kth- 
For CN, we estimate Kth ~ NC, e.g., Kth ~ 12 for N = 4 
and C — 3. Note that k is equal to the dielectric constant 
kq of the medium if the nanotube is immersed in it and is 
equal to (k + l)/2 if the medium is used as a substrate. 
If k < Kth, our theory can still apply at sufficiently low 
n, e.g., in the Wigner crystal regime. 

One possible application of our results for e(n) is a fine- 
tuning of the operational parameters of carbon nanoelec- 
tronic devices |4j, |5( . On a crude level, such devices are 
tiny capacitors made of CN and control metallic gates. 
Precise knowledge of their capacitance per unit length C 
is desirable for their optimal design and efficiency. The 
quantum and many-body effects influence the measured 
value of C according to the equation (see, e.g., Reflioh 

C- 1 = Co" 1 + (K/e 2 ) x -\ X - l = &eldn*, (18) 

where C _1 - (2/k) \n(2D/R) is the inverse classical (ge- 
ometric) capacitance and x' 1 is the inverse thermody- 
namic density of states (ITDOS). The quantum correc- 
tion due to the ITDOS may be nonnegligible if the dis- 
tance D between between the CN and the gate is small 
or if n is low, so that D ~ a. The measurable signature 
of a finite x _1 would be the n-dependence of C. Re- 
cently, capacitance of CN and their junctions was studied 
in Ref. 2QJ by a 3D DFT. It would be interesting to apply 
our theory to the same structures for comparison. 

The sign of the ITDOS is determined by the convexity 
of e(n) curve. From Fig. ^ we see that at low enough 
electron densities ITDOS becomes negative. This phe- 
nomenon is a generic feature of a strongly correlated elec- 
tron matter [21| . Unlike the case of neutral systems, here 
the negative ITDOS does not imply any thermodynamic 
instability but leads instead to a small overscreening of 



an external electric charge. One possible technique to de- 
tect such an overscreening experimentally is the scanned 
probe imaging of the electrostatic potential along an ul- 
trathin wire (e.g., the CN j^) set on a dielectric sub- 
strate. Above the puddles of the electron liquid induced 
by stray random charges, one would see the potential of 
a "wrong" curvature: higher near the center of the pud- 
dle, lower near its ends. The puddles can be intentionally 
created by additional small gates. 

Finally, from e{n) one can extract the n-dependence of 
the LL parameters that influence charge tunneling and 
low-temperature transport in ID wires. Preliminary re- 
sults and their comparison with other work in the liter- 
ature [2^ have been reported in Ref. A more detailed 
investigation that incorporates the results derived in this 
paper will be presented elsewhere 0] . 
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